module montecarlo
contains
    recursive real function factorial(n) result(fact)
        implicit none
        real, intent(in)    :: n
        if(n == 1 .OR. n == 0) then
            fact = 1
            return
        end if
        fact = n*factorial(n-1)
    end function factorial

	subroutine generateEnergyDistribution(n)
        integer, intent(in)   :: n
        integer :: table(n)
        integer :: state(n)
        real    :: energy, maxenergy, minenergy
        integer :: i, j, k
        real, allocatable:: energytable(:,:), energytable_temp(:,:)
        do i = 1, n
            table(i) = 2**(i-1)
        end do
        state = -1
        allocate(energytable(1,2), energytable_temp(1,2))
        energytable = 0
        energytable_temp = 0
        do j = 1, 2**n
            do i = 1, n
                table(i) = table(i) - 1
                if(table(i) == 0) then
                    table(i) = 2**(i-1)
                    state(i) = -1*state(i)
                end if
            end do
            energy =  calculateEnergy(2.0, 1.0, 1.0, state) !(J, H, u, state)
            do k = 1, size(energytable, dim=1)
                if(energy > energytable(k,1) .AND. k < size(energytable, dim=1)) then
                    continue
                elseif(energy > energytable(k,1) .AND. k == size(energytable, dim=1)) then
!                     if(energy>0) print*, "a"
                    deallocate(energytable_temp)
                    allocate(energytable_temp( size(energytable, dim=1) +1, 2))
                    energytable_temp(1:k,:) = energytable(1:k,:)
                    energytable_temp(k+1,1) = energy; energytable_temp(k+1,2) = 1
                    deallocate(energytable)
                    allocate(energytable(size(energytable_temp, dim=1)+1,2))
                    energytable = energytable_temp
                    exit
                end if
                if(energy < energytable(k,1) .AND. k==1) then
!                     if(energy>0) print*, "energy < energytable(k,1) .AND. k==1"
                    deallocate(energytable_temp)
                    allocate(energytable_temp( size(energytable, dim=1) +1, 2))
                    energytable_temp(1,1) = energy; energytable_temp(1,2) = 1
                    energytable_temp(2:,:) = energytable(1:,:)
                    deallocate(energytable)
                    allocate(energytable(size(energytable_temp, dim=1)+1,2))
                    energytable = energytable_temp
                    exit
                elseif(energy < energytable(k,1) .AND. k/=1) then
!                     if(energy>0) print*, "energy < energytable(k,1) .AND. k/=1"
                    deallocate(energytable_temp)
                    allocate(energytable_temp(size(energytable, dim=1)+1, 2))
                    energytable_temp(1:k-1,:) = energytable(1:k-1,:)
                    energytable_temp(k,1) = energy; energytable_temp(k,2) = 1
                    energytable_temp(k+1:,:) = energytable(k:,:)
                    deallocate(energytable)
                    allocate(energytable(size(energytable_temp, dim=1)+1,2))
                    energytable = energytable_temp
                    exit
                end if
                if(energy == energytable(k,1)) then
!                     if(energy>0) print*, "energy == energytable(k,1)"
                    energytable(k,2) = energytable(k,2) + 1
                    exit
                end if
            end do
        end do
        do i = 1, size(energytable, dim=1)
           print*, energytable(i,1), energytable(i,2)/2**n 
        end do
    end subroutine generateEnergyDistribution

    real function calculateEnergy(Jk, Hk, uk, statelinear)
        real, intent(in)                :: Jk, Hk, uk
        integer, intent(in)             :: statelinear(:)
        integer, allocatable    :: state(:,:)
        integer                 :: n
        real                    :: energy
        integer                 :: i, j
        n = sqrt(real(size(statelinear, dim=1))) + 2 ! 2 is added to add buffer lines to make a toriod
        allocate(state(n,n))
        state = 0
        do j = 2, n-1
            do i = 2, n-1
                state(j, i) = statelinear((j-2)*(n-2) + i-1)
            end do
        end do
        !adding buffer rows and columns at the end points
        state(1,2:n-1) = state(n-1,2:n-1)
        state(n,2:n-1) = state(2,2:n-1)
        state(2:n-1,1) = state(2:n-1,n-1)
        state(2:n-1,n) = state(2:n-1,2)
        energy = 0
        do j = 2, n-1
            do i = 2, n-1
                energy = energy + (-Jk/2.0*(state(j-1, i)+state(j, i+1)+state(j+1, i)+state(j, i-1)) - uk*Hk)*state(j,i)
            end do
        end do
        calculateEnergy = energy
    end function calculateEnergy

end module montecarlo

program test
    use montecarlo
    implicit none
    call generateEnergyDistribution(16)
end program
